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ABSTRACT 


This  paper  presents  the  results  of  empirically  testing  eight  alternative 
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multipliers  for  a multiplicative  congruential  generator  with  modulus  2 -1. 

The  LLRANDOM  number  package  (Learmonth  and  Lewis,  1973)  uses  one  of  the  multi- 
pliers, the  simulation  programming  language  SIMSCRIPT  II  uses  a second  and  the 
remaining  six  are  the  best  of  50  candidate  multipliers  studied  by  Hoaglin  (1976) 
using  the  theoretical  spectral  and  lattice  tests.  The  battery  of  tests  raises 
serious  doubts  about  three  of  the  multipliers,  including  the  one  in  LLRANDOM. 

The  power  of  the  tests  is  demonstrated  by  their  rejection  of  RANDU,  a notably 
poor  random  number  generator.  A comparison  of  the  results  for  the  eight  multi- 
pliers with  the  eight  worst  multipliers  (with  regard  to  2-tuples)  in  Hoaglin 
(1976)  failed  to  show  any  apparent  gross  differences.  Since  examination  of 
performance  on  the  lattice  test  revealed  that  the  16  multipliers  clustered  in 
performance  when  compared  to  RANDU,  one  may  conjecture  that  the  poorer  lattice 
test  performance  for  the  worst  eight  is  too  subtle  for  detection  by  our  empiri- 
cal tests.  Since  this  failure  may  not  be  as  serious  as  the  lattice  test  implies. 


one  may  want  to  revise  upward  the  criteria  for  acceptable  performance  on  the 
lattice  tests  and  similarly  on  the  spectral  test.  In  section  9 an  analysis  of 
algorithms  for  generating  random  numbers  reveals  that  the  choice  of  coding  can 


significantly  affect  execution  time  for  different  multipliers. 
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1.  Introduction 

I 

This  paper  presents  the  results  of  an  empirical  evaluation  of  eight 
suggested  values  for  the  multiplier  A in  the  prime  modulus  multiplicative 
congruential  random  number  generator 

(1)  Zi  - AZi-l  (mod  231-1)  . 

Lewis  et  al . (1969)  suggested  multiplier  I and  it  is  in  common  use  in  the 
LL RAN DOM  random  number  package  (Learmonth  and  Lewis,  1973),  in  APL  (Katzan,  1971 ) , 
in  the  simulation  program  language  SIMPL/1  (IBM, 1972)  and  in  the  IMSL  Library 
(1977).  Payne  et  al . (1969)  suggested  multiplier  II,  and  its  most  common  use 
is  in  the  simulation  programming  language  SIMSCRIPT  II.  The  remaining  six 
suggestions  come  from  a study  of  50  multipliers  by  Hoaglin  (1976)  who  evaluated 
them  as  best  using  the  spectral  test  (Coveyou  and  MacPherson,  1967)  and  the  lattice 
test  (Marsagl ia, 1971 ) . It  is  well  known  (Marsagl ia,  1968)  that  all  linear  con- 
gruential generators  are  flawed.  The  spectral  and  lattice  tests  provide  theo- 
retical methods  of  assessing  the  seriousness  of  the  flaw  in  a multiplier. 

These  tests  showed  that  the  six  multipliers  considered  here  had  less  serious 
flaws  than  the  remaining  44  multipliers  studied. 

In  the  present  paper  we  examine  the  performance  of  these  eight  multipliers  > 

when  subjected  to  a battery  of  tests  designed  to  detect  departures  from  random- 
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ness  in  a sequence  of  numbers.  Let  = Z../(2  -1)  for  i = l,...,n. 

Hypotheses  to  be  tested  include: 

Hq.  {Uj}  is  a sequence  of  i.i.d.  random  variables. 

Hj . Uj  has  a uniform  distribution  on  (0,1). 

H2'  ^U2i  i’u2i)  has  a un,ifon,n  distribution  on  the  unit  square. 

Hq.  ^U3i  2,U3i  i »u3i ) has  3 uniform  distribution  on  the  unit  cube. 

H4.  Hq,  H^ , H2  and  Hq  hold  simultaneously.  (If  HQ  and  H1  hold,  H^ 


and  Hq  are  redundant.) 


i 


2 > 

{ 

For  each  multiplier,  the  results  include  a separate  evaluation  of  the  test  of  , 

each  hypothesis  HQ,  H1 , H 2 and  H3  and  a collective  evaluation  via  of  the 

performance  on  the  four  hypotheses  taken  together.  The  collective  evaluation 

reveals  that  three  of  the  candidate  multipliers  (I,  III  and  VII)  are  seriously 

suspect. 

In  any  empirical  study  of  this  sort,  one  hopes  that  the  collective  power 
of  the  tests  employed  is  sufficiently  high  to  detect  departures  from  randomness 
when  they  exist.  We  demonstrate  this  ability  to  detect  departures  in  an  appli- 
cation of  these  same  tests  to  samples  generated  by  RANDU,  the  random  number 
generator  used  until  recently  in  IBM's  Scientific  Subroutine  Library.  Marsaglia 
(1971)  has  demonstrated  the  notably  poor  performance  of  this  generator  on  the 
lattice  test. 

In  a further  attempt  to  demonstrate  the  power  to  detect  nonrandomness,  we 
applied  the  tests  to  the  eight  multipliers  in  Hoaglin's  study  that  performed 
poorest  on  the  spectral  and  lattice  tests  for  2-tuples.  An  analysis  of  the 
results  of  the  best  and  worst  in  the  Hoaglin  study  showed  no  difference  of 

* 

note,  leading  to  several  conjectures.  One  is  that  the  tests  cannot  consistently 
detect  departures  of  the  magnitude  examined.  (However,  they  can  detect  departures 
as  substantial  as  RANDU's.)  Another  is  that  the  criteria  suggested  in  Knuth  (1969) 
and  Marsaglia  (1971)  for  a multiplier  to  "pass"  the  spectral  and  lattice  tests, 
respectively,  may  be  considerably  more  conservative  than  first  thought.  We 
discuss  these  conjectures  in  more  detail  in  Section  8. 

2.  Distribution  Testing  Procedures  ' 

% 

For  each  multiplier,  we  collected  100  independent  samples,  each  with 

n = 200,000  observations.  For  sample  i and  hypothesis  j a statistic  T.  . , 

■ »J 

was  computed.  Then  for  hypothesis  j,  j ,T2  j*****T100  j were  subJected  t0 


3 


a battery  of  tests.  Let  T.  ■ have  cumulative  distribution  function  (c.d.f.) 

■ 

G.  under  hypothesis  j.  Then  P . = 1-G.(T.  .)  has  the  uniform  distribution 

J I »J  J I >J 

on  (0,1)  and 


(2) 


Fn,j(t)  ' n ^i=l  ^O.t]  ^ P i , j ^ 


0 s t < 1 , 


where  I denotes  the  indicator  function,  is  an  empirical  c.d.f.  If  hypothesis 
j is  true 


(3) 


Dn,j  = sJplFn,j(t)  - ‘I 


has  the  Kolmogorov-Smirnov  (K-S)  distribution.  Also 

<4>  V,j  " 

has  the  uniform  distribution  on  (0,1)  (Dwas^l958).  Finally,  for  large  n 

(5)  = n J1  ftFn  j (t)-tJ2/t(l-t)}dt 


has  a distribution  given  by  Anderson  and  Darling  (1952). 

2 

The  motivation  for  goodness-of-fit  tests  based  on  D„  , V„  . and  A„  . 

3 n,j  n,j  n,j 

arises  from  the  distinct  departures  from  behavior  under  hypothesis  j that 

each  is  designed  to  detect.  The  quantity  D , measures  the  maximal  absolute 

n»j 

deviation  between  the  empirical  and  the  hypothesized  c.d.f. ‘s,  V . measures 

n » j 

2 

the  proportion  of  F . that  lies  below  the  hypothesized  c.d.f.  and  A . 

n , j n , j 

measures  the  extent  of  deviation,  principally  in  the  tails  of  the  empirical 
c.d.f. 


3.  Testing  for  Independence  HQ 

To  test  Hq  we  relied  on  runs-up-and-down  statistics.  Let  Ri  k and 
k denote  the  numbers  of  runs  up  and  down,  respectively,  of  length  k on 


I 


replication  i.  Then  the  statistics  for  HQ  are 

<6>  Ti,<r  ^t=1MRi,k  - E(Ri,kH[Ri,<  - 

*dMtRi,k-E'Ri,k)][Si.«  ~e<rm» 

+ cMtRitk  * E<si,k>i[RM  - E<Ri,*»t  < ■ '•••••">» 

and  each  asymptotically  has  the  chi-squared  distribution  with  12  degrees  of  free- 
dom. The  quantities  c^  and  d^  are  computable  from  the  covariance  matrix 
of  Ri  i , . . . , Ri  1,...,Ri  g (Levene,  1953) . For  200,000  observations  one 

has 


(7) 


e ( R-j 

,1>  = 

ic£~ 

LU 

.i>- 

41670 

E (Ri 

,2}  = 

E(Ri 

,2}  = 

18330 

E(Ri 

,3)  = 

E(Ri 

,3)  = 

5280 

E(Ri 

,4}  = 

E(Ri 

,4}  = 

1150 

E(Ri 

,5}  = 

E(^i 

,5}  = 

200 

E(Ri 

,6}  = 

E(*1 

,6}  = 

30 

E(Ri 

,7+)  = 

= E(R 

i,7+) 

= 4 

Here  R.  7+  and  R.  7+  denote  the  numbers  of  runs  up  and  down,  respectively, 
of  length  7 or  more.  Although  a statistic  similar  to  (6)  that  incorporates 
either  R^  or  can  be  constructed  and  asymptotically  has  a chi -squared 

distribution  with  13  degrees  of  freedom,  the  small  value  of  E(R.  7+)  in  the 
present  case  encouraged  us  to  work  with  (6)  to  avoid  any  discretization  error 
that,  say,  R^  j+  might  induce.  Note  that  inclusion  of  R^  and  R^  7+ 
in  (6)  would  produce  a degenerate  distribution  for  T..  Q (Wolfowitz,  1944). 

By  way  of  interpretation,  excessive  short  (long)  runs  imply  more  mixing 


It 


5 


(clustering)  than  one  would  expect  to  find  in  a purely  random  sequence.  The 
decision  to  examine  runs  up  and  runs  down  explicitly  rather  than  study  their 
sum  was  motivated  by  the  finding  in  Tootill  et  al . (1971)  that  runs  up  and 
runs  down  can  exhibit  distinct  behavioral  patterns  for  certain  Tausworthe 
random  number  generators  (Tausworthe,  1965) . Although  the  generator  (1)  is 
not  of  this  type,  we  decided  to  allow  for  the  possibility  of  analogous  behavioral 
distinctions. 

Column  3 of  Table  1 lists  the  P-values  for  the  test  statistics  (3),  (4) 

and  (5)  for  each  multiplier.  In  particular,  a P-value  for  the  K-S  (A-D)  test 

is  the  probability  that  a random  variable  from  the  K-S  (A-D)  distribution 
2 

exceeds  D10Q  Q (A-jqq  q)  in  value  under  HQ.  For  (4)  the  P-value  is 
2 min(V100  1-1^  Q0  q),  since  V]oo  j is  a uniform  deviate  under  H ^ . Except 
for  RANDU,  the  P-values  in  column  3 give  little  cause  for  concern. 


4.  Testing  for  Uniformity 

To  test  Hj  we  chose  a chi-squared  goodness-of-fit  statistic.  Consider 
K cells  on  the  unit  interval  each  of  length  1/K.  Let  ^ denote  the  number 
of  the  N observations  on  replication  i that  fall  into  the  interval 
((k-l)/K,  k/K].  Then  under 


(8) 


£,  <Ni,k-N/K>2  = £ t,  "U-N 


i = 1 100 


each  asymptotically  have  a chi-squared  distribution  with  K-l  degrees  of  freedom. 
Choosing  K = 2^  = 4096  implied  a cell  width  of  1/K  = 0.000244140625  and 
enabled  us  to  test  the  first  12  bits  of  U.  . 

2 

Column  4 of  Table  1 shows  the  P-values  for  D^qq  ^ , V^QOj  and  A^QQ  1 
for  each  multiplier.  The  P-values  again  arouse  little  suspicion.  The  failure 
to  detect  nonuniformity  in  RANDU  may  be  an  indication  of  the  poor  power  of  this 


P-Values  for  Testing  Hypotheses 
(RANDU  and  Multipliers  I through  VIII) 


I 

6 , 


P-Values  for  Testing  Hypotheses 
(RANDU  and  Multipliers  I through  VIII) 
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test  or  an  indication  that  RANDU  has  the  local  property  of  one-dimensional 
uniformity. 


1 


5.  Testing  for  Bivariate  Uniformity  H2 


Hypothesis  H2  is  designed  to  detect  nonuniformity  when  the  U.  are 
taken  in  nonoverlapping  pairs  or  2-tuples.  One  motivation  for  this  testing 
arises  from  the  theoretical  observation  in  Marsaglia  (1968)  that  the  random- 
ness of  k-tuples  becomes  more  suspect  as  k increases.  The  spectral  and 
lattice  tests  support  this  observation.  In  particular,  see  Hoaglin  (1976)  and 
Marsaglia  (1971).  Ideally  one  would  like  to  test  for  the  uniformity  in  distri- 
bution of  k-tuples  of  the  k-dimensional  unit  hypercube.  In  practice,  such 
testing  is  excessively  expensive,  even  for  k = 2. 

Let  us  divide  the  unit  interval  into  K cells,  each  of  width  1/K.  Let 
denote  the  frequency  with  which  nonoverlapping  2-tuples  fall  into  the 
square  (((j-l)/K,  j/K]  *((k-l)/K,  k/K] ) on  replication  i.  For  fixed  K the 
quantities 


(9) 


T i » 2 


1C 


= y 

N L- 


j,k=l 


(Njjk  - N/K2)2 


1 100 


2 

each  asymptotically  have  the  chi-squared  distribution  with  K - 1 degrees  of 

freedom.  Suppose  we  had  chosen  as  before  K = 4096.  Then  there  would  be 
2 

K = 16777216  cells.  To  guarantee  a mean  of  5 per  cell  under  H2  would 

require  N > 80  million  observations  per  replication  or  over  8 billion  observa- 

31 

tions  per  multiplier.  Since  2 -1  < 4.3  billion,  such  a sample  size  is  not 

possible  using  this  test  procedure.  Because  of  this  demonstrated  excessiveness, 

we  chose  K = 128,  which  required  16384  cells,  implied  a cell  width  of  0.0078125 

2 

and,  for  N = 200,000,  a mean  of  n/K  = 12.21  per  cell  under  W^.  This  choice 
of  K enabled  us  to  study  the  first  7 bits  of  each  coordinate  of  a 2-tuple. 


9 


Column  5 in  Table  1 disp’ays  the  P-values  for  D1Q0  2,  V1Qp  2 and 
A100  2 ‘ ^ though  V1Q0  2 appears  low  for  multipliers  I and  III,  we  do  not 
regard  this  suspiciously  at  this  point,  since  we  are  simultaneously  examining 

2 

3x8  = 24  P-values.  What  is  of  concern  to  us  is  when  D-jqq  j,  V^q  j and  A^g 
all  give  low  P-values  for  a given  multiplier.  We  demonstrate  cases  of  this  sort 
shortly  when  testing  H^.  Note  the  extremely  poor  showing  for  RANDU. 


6.  Testing  for  Trivariate  Uniformity 

Hypothesis  tests  for  nonuniformity  when  the  U^  are  taken  in  nonover- 
lappinq  triplets  or  3-tuples.  The  motivation  is  similar  to  that  for  testing 
H2-  Let  denote  the  frequency  on  replication  i with  which  nonoverlap- 

ping triplets  fall  into  the  cube  ( ( (j-1 )/K,j/K]  x((k-l)/K,  k/K]  x ( (£-1 )/K,£/K]) 


for  j,k,£  = 1 , . . . ,K  with  K specified.  Then 

<,0)  Ti,i'T  !lu=,'V-w3)Z 


i = 1 100 


each  asymptotically  have  the  chi-squared  distribution  with  K -1  degrees  of 

4 3 

freedom.  Here  we  choose  K = 2 =16  which  gives  K = 16x16x16  = 4096  cells 

3 

and  an  expectation  of  N/K  = 48.83  per  cell.  A perusal  of  the  P-values  in 
column  6 of  Table  1 indicates  a rejection  of  Hg  for  RANDU  and  reveals  some 
concern  about  multiplier  III  and  a lesser  concern  about  multiplier  II. 


7.  An  Omnibus  Test  H4 

In  reviewing  the  P-values  for  HQ,  H^ , H2  and  Hg  one  finds  it  difficult 
to  say  that  any  multipl ier except  RANDU  seriously  errs  in  an  omnibus  sense.  This 
finding  is  entirely  in  keeping  with  other  empirical  investigations  using  tests  with 


vague  alternative  hypotheses.  To  overcome  this  inadequacy,  we  test  H^:  H^, 

H.| , and  H^  are  true  simultaneously. 


i 
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Recall  that  T.  n,  T.  .,  T.  9,  T.  7 are  the  statistics  on  replication  i 

1 >U  1 ) I 1 1 

for  hypotheses  0 through  3.  Under  H,  P,  , = 1 -G  - (T . .)  for  j = 0,...,3 

^ I J IjJ 

each  have  the  uniform  distribution  on  (0,1).  Let 


zi,j  * 


-1 


where  (•)  denotes  the  inverse  function  of  the  unit  normal  c.d.f.  Let 


01) 


Zi,min  " min(Zi,0’  Zi,l*  Z i , 2 ’ Zi,3) 


Zi,max  = max(Zi,0’  Zi .1  * Zi,2’  Zi,3) 


Under  the  statistics 


(12) 


Ti,4  l-'K-Z..  ,min’  -Zi,min’  ~Zi,min’  Zi,min^ 


★ 

Ti,4  ” ^-<I>^Zi .max’  Zi,max’  Zi,max’  Zi,max  ’ 


where  <t(*, ’,',•)  denotes  a four-dimensional  multivariate  normal  distribution 
each  have  the  uniform  distribution  on  (0,1).  Under  H^,  T\  ^ is  the  probabil- 
ity of  observing  a minimum  P-value  smaller  than  min(Pi  g,  P..  -j,  P^  2»  p-j  3)  and 
★ 

T.  . is  the  probability  of  observing  a maximum  P-value  greater  than 
max(Pi,0’  Pi,r  Pi ,2’  Pi,3)* 

Figure  1 

Possible  Arrangements  of  Pi  q’  Pi  1*  Pi  2 and  Pi  3 


L 


J 


(a) 


44  •-  » 


(b) 
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★ 

Figure  1 provides  a rationale  for  studying  T.  . and  T.  . . Figure  la 

shows  three  "reasonable"  P.  . 's  and  one  excessively  small  one.  One  way  to 

' »J 

evaluate  the  extent  to  which  this  behavior  is  consistent  with  H4  is  to  test 

Tl  4 ^or  un''formity  on  (0*1)-  Figure  lb  shows  a case  in  which 

no  P.  n,  P.  ,,  P.  9,  or  P.  is  small  nor  is  any  moderately  large.  Here  one 
★ ★ 

tests  T1  4 ^or  uniformity  on  (0.1)  to  determine  the  extent  to 

which  the  behavior  in  Figure  lb  is  consistent  with  H^.  Although  tests  that 
explicitly  include  the  other  order  statistics  can  be  performed,  we  felt  that 
the  two  selected  tests  represented  sufficiently  extreme  situations  that,  if  sig- 
nificant, would  be  regarded  as  unequivocal  evidence  for  disqualifying  a multiplier. 

Because  P.  n,  P.  ,,  P.  9 and  P.  , are  based  on  the  same  sample  data  of 
200,000  observations,  treating  them  as  independent  would  be  unreasonable.  Instead 
we  concentrate  on  Q,  Z.  Z^  2,  and  Z.  3 to  acknowledge  the  dependence 
appropriately.  Under  these  quantities  have  a four-dimensional  multivariate 
normal  distribution  with  mean  vector  0 and  unknown  covariance  matrix  e_  . Con- 
sider the  Cholesky  decomposition  representation 


Then  one  can  express  (12)  as 


12 


Since  I is  unknown,  the  computation  ot  (14)  is  not  directly  possible.  Here 
we  estimated  £ by 


°ij  n-1  ^k=i  (Zk,i"Zi^Zk,j"V 


Zi  = n l 


kal  k.i 


n = 200,000 


from  which  we  estimated  Cg,  c^ , and  c-j  and  (14)  by  substituting 

★ 

the  estimates  for  Cg,  c^ , and  c^.  If  Ti  4 and  T..  ^ have  the 

proper  distribution  under  then  estimating  c will  cause  the  data  to  appear 

to  follow  more  closely  than  using  the  true  c . Thus  the  tests  employed 
here  are  conservacive. 

Columns  7 and  8 of  Table  1 show  the  P-values  for  H^.  Here  the  test 
based  on  the  maximum  provides  strong  evidence  against  multiplier  I and  the 
test  based  on  the  minimum  provides  strong  evidence  against  multipliers  III  and 
VII.  RANDU's  failure  on  the  minimum  test  is  consistent  with  its  poor  showing  on 
Hg,  H^  and  Hj. 

In  interpreting  the  results  for  the  minimum  or  maximum  test  in  column  7 
or  8,  one  needs  to  pay  attention  to  the  issue  of  multiplicity  across  the  eight 
multipliers  (ignoring  RANDU).  Under  the  hypothesis  that  the  eight  multipliers 
each  satisfy  H^,  the  probability  of  observing  at  least  two  K-S  test  statistics 
there  with  P-values  £ .006  on  the  minimum  test  is  .0010.  Alternatively,  under 
the  same  hypothesis,  the  probability  of  observing  at  least  one  K-S  test  statistic 

O 

with  P-value  s .004  on  the  maximum  test  is  1 - ( . 996)  = .0316.  We  regard  these 
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probabilities  as  sufficiently  small  to  reject  equality  of  behavior  on  the  minimum 
and  maximum  tests  across  the  eight  multipliers.  Therefore,  we  regard  the  poor 
behavior  of  multipliers  I,  III  and  VII  as  real  and  not  simply  due  to  sampling. 

8.  Other  Prime  Modulus  Multipliers 

As  mentioned  earlier,  the  six  multipliers  selected  for  study  from  Hoaglin's 
paper  were  evaluated  there  as  the  best  performers  on  the  spectral  and  lattice 
tests.  A question  naturally  arises  as  to  how  sensitive  the  battery  of  empirical 
tests  presented  here  is  to  the  flaws  present  in  the  remaining  44  prime  modulus 
multipliers  there.  As  a first  inquiry,  we  selected  the  worst  eight  multipliers 
there,  based  on  their  performance  on  the  spectral  and  lattice  tests  for  2-tuples, 
and  subjected  them  to  the  same  battery  of  empirical  tests.  The  results  in  Table 
2 show  no  apparent  gross  departure  from  those  in  Table  1. 

To  understand  why  the  tests  detect  the  poor  behavior  of  RANDU  but  do  not 
, reject  all  multipliers  IX  through  XVI,  we  return  to  lattice  test  performance. 

\ 

The  lattice  test  relies  on  the  reduced  basis  vectors  of  the  lattice  of  n-tuples. 

This  is  the  basis  that  is  as  close  as  possible  to  being  orthogonal,  given  the 
constraints  of  the  lattice  (Hoagl in, 1976) . The  ratio  Ln  of  the  lengths  of  the 
longest  to  the  shortest  basis  vectors  provides  the  figure  of  merit.  A ratio  of 
unity  is  ideal. 

Table  3 lists  this  ratio  for  all  multipliers  examined  in  this  study.  As 
measured  by  Ln,  RANDU  performs  considerably  more  poorly  for  n = 3,  4 and  5. 

Therefore,  one  possible  explanation  of  the  less  than  complete  rejection  of  mul- 
tipliers IX  through  XVI  is  that  the  flaws  revealed  by  the  lattice  test  are 
t too  subtle  for  detection.  Since  the  rule  of  thumb  for  considering  a multiplier 

acceptable  with  regard  to  n-tuples  is  Ln  < 2,  one  may  want  to  revise  this 
bound  upward  in  view  of  our  results.  A similar  upward  revision  may  apply  to 
the  acceptance  criteria  for  the  spectral  test. 

1 


P-values  for  Testing  Hypotheses 
(Multipliers  IX  through  XIV) 
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Table  3 

Lattice  Test  Performance3 


Multipl ier 

L2 

L3 

L4 

L5 

L6 

RANDU 

1.00 

1819 

1872 

274 

— 

I 

7.60 

3.39 

2.07 

1.67 

3.36 

II 

1.26 

2.92 

1.64 

1 .52 

2.11 

III 

1.03 

2.11 

1 .88 

1.55 

1.13 

IV 

1 .07 

1.34 

1.17 

1.62 

1.89 

V 

2.82 

2.63 

1.50 

1 .20 

1.86 

VI 

1.88 

2.84 

1.39 

1.51 

1.24 

VII 

1.64 

2.36 

1.71 

1.24 

1.24 

VIII 

1 .69 

1.62 

1.88 

1.38 

1.41 

IX 

27.91 

1.20 

1.71 

2.00 

1.21 

X 

21.94 

3.10 

1.45 

3.78 

2.47 

XI 

18.1 

1.17 

3.14 

1.41 

2.37 

XII 

12.16 

3.48 

1.54 

1.41 

1.33 

XIII 

10.81 

4.30 

2.84 

2.73 

1.32 

XIV 

8.77 

4.02 

1.28 

2.65 

1.30 

XV 

8.33 

1.19 

2.32 

2.10 

2.09 

XVI 

7.05 

2.53 

1.95 

3.18 

1.59 

3 Source:  For  RANDU,  (Marsaglia  1971).  For  Multipliers  I through  XVI,  (Hoaglin  1976). 
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9.  Execution  Time 

In  Table  1 we  listed  the  relative  mean  execution  times  observed  for  each  mui-  . 

tiplier  (except  RANDU).  These  samples  were  obtained  by  using  the  LLRAi'JDOH  package 
(Learmonth  and  LewiSjl973)  with  a substitution  of  multipliers  for  16807.  The 
apparent  increase  in  execution  time  as  A increases  in  magnitude  makes  the 
desirability  of  large  multipliers  questionable.  To  put  this  issue  in  perspec- 
tive, we  examine  an  algorithm  due  to  Payne  et  al . (1969)  that  is  in  common  use 
for  generating  deviates  from  prime  modulus  generators. 

Suppose  one  has  storage  areas  A,  X and  Y each  capable  of  holding  p 
binary  digits  and  a working  area  D capable  of  holding  2p  binary  digits. 

Payne  et  al . (1969)  developed  the  following  algorithm  with  A containing  the 
multiplier  and  X containing  the  most  recently  generated  random  number. 

LCG1 

1.  D ■*-  X • A . 

2.  Y •*-  i.p.[D/2p_1]  . 

3.  X D - Y • 2P_1  . 

4.  X X + Y . 

5.  Y i.p.[X/2p-1]  . 

6.  If  Y * 0 , X +-  X - ( 2P_  1 -1 ) . 

7.  Return  X . 

Steps  2 and  3 may  be  performed  using  shift  operations  on  D.  Step  5 and 
the  check  in  step  6 may  be  performed  on  p bit  word  machines  by  an  overflow 
check. 

It  is  instructive  to  examine  how  often  the  check  in  step  6 is  satisfied. 

Of  the  possible  (2p~^-2)  different  initial  values  of  X,  A of  these  cause 
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the  check  in  step  6 to  be  satisfied  where 

* i 

i I 

i . p . [ ( A-2 ) (A-l )/2A]  if  A < 2P'2 

i.p.  j(A-2)  (A-l  )/2A  + 2-  A]  if  A > 2P"2  • 

Therefore,  the  average  execution  time  for  LCG1  tends  to  follow  the  relationship 

e0  + e1A*/(2p_1-2) 

where  Bg  is  the  time  to  perform  steps  1 through  7 of  LCG1  without  the 

check  in  step  6 being  satisfied  and  (Sj  is  the  time  to  perform  the  operation 

★ 

in  step  6 when  the  check  is  satisfied.  For  most  values  of  A,  A is  very 
close  to  A/2  and  the  average  execution  time  for  LCG1  is  Bg  + B-|A/2P  . 

In  practice,  execution  time  is  code-dependent  as  well  as  algorithm-dependent. 

In  particular,  a code  that  produces  a very  small  B-j  is  most  desirable.  The 
code  provided  in  Payne  et  al . (1969)  and  used  with  slight  modification  in  SUB- 
SCRIPT II  should  produce  a very  small  b-j  . The  code  is  self-contained.  How- 
ever, the  code  in  LLRANDOM  for  step  6 will  produce  a relatively  large  value 
for  6-|  • This  is  due  to  calls  to  a routine  outside  the  generator  code  (when 
an  overflow  occurs).  Regrettably,  this  observation  did  not  come  to  our  attention 
until  all  runs  were  completed.  In  summary,  all  implementations  of  the  Payne  et  al . 
algorithm  have  execution  times  proportional  to  A . However,  the  choice  of  code 
determines  the  constant  of  proportionality. 
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